El modelo de regresión lineal múltiple generaliza el modelo lineal simple por incluir \(k\) variables regresoras. Se puede formular de dos formas distintas:
\[ Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_{p-1} X_{ip-1} + \varepsilon_i,\quad \varepsilon_i \sim N(0, \sigma^2) \label{eq:modelo} \]
donde:
\(Y_i\) es la variable respuesta para la observación \(i\), con \(i=1,\ldots,n\).
\(\beta_j\) son los coeficientes del modelo, con \(j=0,\ldots,p-1\).
\(X_{ij}\) es la \(j\)-ésima covariable para esa observación.
\(\varepsilon_i\) es el término de error aleatorio.
\(p\) es la cantidad de parámetros a estimar incluyendo el intercept.
En forma matricial, el modelo se escribe como:
\[ \mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \]
donde:
\(\mathbf{Y}\) es un vector \(n \times 1\) de respuestas,
\(\mathbf{X}\) es una matriz \(n \times p\) de diseño (incluye una columna de 1s),
\(\boldsymbol{\beta}\) es un vector \(p \times 1\) de parámetros,
\(\boldsymbol{\varepsilon}\) es un vector de errores aleatorios, con \(\boldsymbol{\varepsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I})\).
Estos vectores y matrices son:
\[\begin{equation} \begin{array}{cc} \underset{n \times 1}{\mathbf{Y}}=\left[\begin{array}{c} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{array}\right] & \underset{n \times p}{\mathbf{X}}=\left[\begin{array}{ccccc} 1 & X_{11} & X_{12} & \cdots & X_{1, p-1} \\ 1 & X_{21} & X_{22} & \cdots & X_{2, p-1} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & X_{n 1} & X_{n 2} & \cdots & X_{n, p-1} \end{array}\right] \\ \underset{p \times 1}{\boldsymbol{\beta}}=\left[\begin{array}{c} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_{p-1} \end{array}\right] & \underset{n \times 1}{\varepsilon}=\left[\begin{array}{c} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{array}\right] \end{array} \end{equation}\]
Notemos que, para el caso de dos variables regresoras, el modelo resulta:
\[ Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2}+\varepsilon_i,\quad con \quad \varepsilon_i \sim N(0, \sigma^2) \]
En este caso, la esperanza condicional de \(Y\) dada las variables \(X_1\) y \(X_2\) es
\[\begin{equation} E(Y|X_{1},X_{2}) = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2}, \label{eq:modelo2} \end{equation}\]
A esta ecuación se la llama superficie de respuesta y su gráfico es un plano. En el caso de regresión lineal simple tenemos que \(E(Y|X) = \beta_0 + \beta_1 X\) y su gráfico es una recta.
\(\operatorname{los} \varepsilon_i\) tienen media cero, \(E\left(\varepsilon_i\right)=0\).
\(\operatorname{los} \varepsilon_i\) tienen todos la misma varianza desconocida que llamaremos \(\sigma^2\) y que es el otro parámetro del modelo, \(\operatorname{Var}\left(\varepsilon_i\right)=\sigma^2\).
\(\operatorname{los} \varepsilon_i\) tienen distribución normal.
\(\operatorname{los} \varepsilon_i\) son independientes entre sí, e independientes de las covariables \(X_{i 1}, X_{i 2}\), \(\ldots, X_{i p-1}\).
\[ \hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y} \]
\[ \min_{\boldsymbol{\beta}} \| \mathbf{y} - \mathbf{X} \boldsymbol{\beta} \|^2 \]
De la misma manera que sucede en regresión lineal simple, se puede probar que el estimador de máxima verosimilitud para los parámetros del modelo coinciden con los de mínimos cuadrados, siempre que los errores se distribuyan en forma normal y sean independientes.
\[ \hat{\sigma}^2=\dfrac{\sum_{i=1}^n (Y_{i}-\widehat{Y}_{i})^2}{n-p}\] donde \(p\) es la cantidad de parámetros a estimar.
Un embotellador de bebidas gaseosas quiere mejorar la eficiencia de sus rutas de distribución. Para eso, está estudiando cuánto tiempo le toma a un representante de ruta atender las máquinas expendedoras en una tienda. Este proceso de atención incluye reponer productos embotellados en la máquina y realizar tareas básicas de limpieza y mantenimiento. Un ingeniero industrial que colabora en el estudio propone que el tiempo que se tarda en atender una máquina depende principalmente de dos factores: la cantidad de cajas de producto que se reponen (\(X_1\)) y la distancia que el representante tiene que caminar dentro de la tienda (\(X_2\) medida en metros). Se han recogido 25 observaciones sobre el tiempo de servicio.
Este ejemplo se obtuvo del libro Introducción al análisis de regresión lineal de Montgmoery et al. que se encuentra en las referencias.
library(ggplot2)
library(GGally)
ggpairs(Tiempos,
columns = 1:3,
mapping = aes(alpha = 0.9),
upper = list(continuous = wrap("points", color = "darkred")),
lower = list(continuous = wrap("points", color = "blue")),
diag = list(continuous = wrap("densityDiag", fill = "lightblue")))scatterplot3d(x=Tiempos$cajas, y=Tiempos$distancia, z=Tiempos$tiempo, pch=16, cex.lab=1,
highlight.3d=TRUE, type="h", xlab="Cantidad de cajas",
ylab="Distancia", zlab="Tiempo")Se puede observar que, a medida que la cantidad de cajas y la distancia aumentan, los tiempos también aumentan.
Se puede realizar un diagrama de dispersión con movimiento.
\[ Y_i=\beta_0+\beta_1 X_{1i}+\beta_2X_{2i}+\varepsilon_i,\quad con \quad \varepsilon_i \sim N(0, \sigma^2) \] donde
\(Y\): tiempo de atención a las máquinas expendedoras.
\(X_1\): cantidad de cajas entregadas.
\(X_2\): distancia recorrida por el representante.
##
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9508 -1.7810 0.0006 1.8293 4.3545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.137487 0.986207 5.209 3.66e-05 ***
## cajas 1.526082 0.134819 11.320 2.12e-10 ***
## distancia 0.008508 0.002955 2.879 0.00897 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.516 on 21 degrees of freedom
## Multiple R-squared: 0.943, Adjusted R-squared: 0.9376
## F-statistic: 173.7 on 2 and 21 DF, p-value: 8.666e-14
\[ E(Y|X_{1},X_{2}) =2.341231+1.615907X_1+0.014385 X_2 \]
# Se crea el grafico 3d y se guarda en un objeto, por ejemplo mi_3d
mi_3d <- scatterplot3d(x=Tiempos$cajas, y=Tiempos$distancia, z=Tiempos$tiempo, pch=16, cex.lab=1,highlight.3d=TRUE, type="h", xlab="Cantidad de cajas",ylab="Distancia ", zlab="Tiempo")
# Para agregar el plano usamos $plane3d( ) con argumento modelo ajustado
mi_3d$plane3d(modelo, lty.box="solid", col="mediumblue")\(\beta_0=2.341231\) indica el tiempo promedio estimado de servicio cuando:
no se entregan cajas \(X_1=0\),
no se camina distancia alguna \(X_2=0\).
\(\beta_1 = 1.615907\) indica que el tiempo de servicio aumenta, en promedio, aproximadamente \(2.34\) minutos cuando se agrega una caja adicional y se mantiene la distancia constante.
\(\beta_2 = 0.014385\) indica que el tiempo de servicio aumenta, en promedio, aproximadamente \(0.01\) minutos cuando la distancia recorrida aumenta un metro y se mantiene la cantidad de cajas constante.
ggplot(data = Tiempos, aes(x = fitted(modelo), y = rstandard(modelo))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "valores ajustados", y = "residuos estandarizados",
title = "Residuos estandarizados vs Valores predichos")Estos gráficos sirven para detectar la existencia de una relación curvilínea entre los residuos y cada variable regresora. Esto indicaría que debiera incluirse un término \(X_i^2\) o \(ln(X_i)\) o alguna otra transformación.
ggplot(data = Tiempos, aes(x = cajas, y = rstandard(modelo))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "cantidad de cajas", y = "residuos estandarizados",
title = "Residuos estandarizados vs Cantidad de cajas")ggplot(data = Tiempos, aes(x = distancia, y = rstandard(modelo))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "distancia", y = "residuos estandarizados",
title = "Residuos estandarizados vs Distancia")residuos.df=data.frame(resid=rstandard(modelo))
ggplot(data = residuos.df, aes(sample = rstandard(modelo))) +
stat_qq() +
stat_qq_line(color = "red") +
labs(title = "QQ Plot de los residuos",
x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
theme_minimal()ggplot(residuos.df, aes(y = rstandard(modelo))) +
geom_boxplot(fill = "skyblue") +
labs(title = "Boxplot de residuos estandarizados",
y = "Residuos estandarizados") +
theme_minimal()\[ R^2=1-\frac{\text { Suma de cuadrados del residuo ($SSRes$) }}{\text { Suma total de cuadrados (SST) }} \] - Este coeficiente tiene la misma interpretación que en regresión lineal simple.
Pero…, hay que notar que si agregamos variables al modelo, el valor de \(SSRes\) disminuye.
Es decir, si agregamos variables al modelo, el valor de \(R^2\) aumenta en forma espúrea.
Por este motivo,cuando queremos comparar modelos de regresión que tienen distinta cantidad de variables, no conviene usar el \(R^2\),
Conviene usar otra medida que corrige ese valor teniendo en cuenta la cantidad de variables que se usaron en el modelo.
El \(R^2_a\) penaliza por el número de variables en el modelo y solo aumenta si la nueva variable mejora realmente el modelo. Se define:
\[ R_{a}^2=1-\left(\frac{\left(1-R^2\right)(n-1)}{n-k-1}\right) \]
donde:
\(n\) es el número de observaciones,
\(k\) es el número de predictores.
Se puede analizar si cada parámetro es indivdualmente significativamente distinto de cero, de la misma forma que hicimos en regresión lineal simple. Los test \(\textit{t}\) sirven para contrastar si cada parámetro individual \(\textit{\beta_j}\) es significativamente distinto de 0. Es decir, testear las siguientes hipótesis:
\[ H_0: \ \beta_j=0 \text{ vs. } \ H_1: \ \beta_j \neq 0 \ \text{ con } \ j=1,2\] Donde:
\(\beta_1\) es el efecto de la variable \(\textit{cajas}\) (que mide el número de productos entregados).
\(\beta_2\) es el efecto de la variable \(\textit{distancia}\) (que mide la distancia caminada).
\[T=\dfrac{\hat{\beta}_j}{se(\hat{\beta}_j)} \sim t_{n-p-1} \ \text{ bajo } \ H_0\]
Se rechaza \(H_0\) si \(|t_c|>t_{\alpha/2,n-p-1}.\)
Es importante notar que esta es una prueba parcial o condicional, ya que el coeficiente de regresión \(\beta_j\) depende de la presencia de todas las demás variables explicativas \(x_i\) (con \(i \neq j\)) en el modelo. Por lo tanto, el test \(t\) para \(\beta_j\) evalúa el aporte de la variable \(x_j\) al modelo, condicional a que las demás variables ya están incluidas. En otras palabras, nos dice si \(x_j\) mejora significativamente la explicación de la variable respuesta, dado que las demás regresoras ya están en el modelo.
##
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9508 -1.7810 0.0006 1.8293 4.3545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.137487 0.986207 5.209 3.66e-05 ***
## cajas 1.526082 0.134819 11.320 2.12e-10 ***
## distancia 0.008508 0.002955 2.879 0.00897 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.516 on 21 degrees of freedom
## Multiple R-squared: 0.943, Adjusted R-squared: 0.9376
## F-statistic: 173.7 on 2 and 21 DF, p-value: 8.666e-14
Todos los p-valores son chicos, quiere decir que se rechaza la hipótesis nula a cualquier nivel indicando que cada variable aporta información explicativa más allá del azar.
Por lo tanto, tanto la variable \(\textit{cajas}\) como la variable \(\textit{distancia}\) influyen significativamente en el tiempo de entrega.
Si algún p-valor es grande, no se rechaza la hipótesis nula indicando que la/s variable/s no influyen significativamente en el modelo, pudiendo también indicar presencia de colinealidad (ver VIF más abajo).
Este test permite evaluar si el modelo es globalmente significativo, es decir, si al menos una de las variables regresoras explica la variabilidad de la respuesta.
Este test se basa en el cociente de la mejora debida al modelo (SSReg) y la diferencia entre el modelo y los datos observados (SSRes).
El objetivo es ver si al menos uno de los coeficientes \(\beta_j\) de los predictores es significativamente distinto de cero.
\[H_0: \beta_1=\ldots=\beta_{p}=0 \text{ vs. } H_1: \beta_j=0 \text{ para algún } j=1,\ldots,p\]
\[SS_T=SSReg+SSRes \]
\[\text{F}=\dfrac{\dfrac{\text{Variación explicada por el modelo}}{p-1}}{\dfrac{\text{Variación residual}}{n−p}}=\dfrac{{\dfrac{SSReg}{p-1}}}{{\dfrac{SSRes}{n-p}}}=\dfrac{MSReg}{MSRes}.\]
Donde:
\[\text{F}\sim F_{p-1,n-p}\]
Regla de decisión:
Valores grandes de F corresponde a p-valores pequeños.
Si se rechaza \(H_0\) indica que no todos los coeficientes del modelo son nulos.
Un valor alto de \(F\) indica que el modelo es significativo.
##
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9508 -1.7810 0.0006 1.8293 4.3545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.137487 0.986207 5.209 3.66e-05 ***
## cajas 1.526082 0.134819 11.320 2.12e-10 ***
## distancia 0.008508 0.002955 2.879 0.00897 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.516 on 21 degrees of freedom
## Multiple R-squared: 0.943, Adjusted R-squared: 0.9376
## F-statistic: 173.7 on 2 and 21 DF, p-value: 8.666e-14
En el último renglón del summary vemos que
Por lo tanto el pvalor es muy chico, se rechaza \(H_0\) a nivel \(5\%\) indicando que, al menos algún \(\beta_j\) es no nulo.
Recordemos
\[ p-\text { valor }=P\left(F_{p-1, n-p}>F_{obs}\right) \]
## [1] 24
## [1] 8.648615e-14
al ajustar el modelo estamos usando \(p\) valores (los coeficientes beta) que se determinan a partir de los datos.
esos valores se obtienen resolviendo \(p\) ecuaciones, que se conocen como las ecuaciones normales.
estas ecuaciones surgen al minimizar el error del modelo: se derivan \(p\) veces (una por cada parámetro) y se igualan a cero.
entonces, los residuos no son completamente libres: están “condicionados” por esas \(p\) ecuaciones.
luego, si conocemos \(n-p\) residuos, los otros \(p\) se pueden calcular a partir de ellos.
por eso, decimos que los residuos tienen \(n-p\) grados de libertad.
En un modelo de regresión lineal múltiple, los coeficientes estimados representan el cambio esperado en la variable respuesta asociado con una unidad de cambio en cada predictor, manteniendo los demás constantes. Sin embargo, estas estimaciones son sujetas a variabilidad muestral.
Para cuantificar la precisión de estos estimadores, se calculan intervalos de confianza (IC), que representan rangos plausibles donde se espera que se encuentre el valor real del parámetro con un cierto nivel de confianza (usualmente 95%).
Matemáticamente, un intervalo de confianza al 95% para un coeficiente \(\beta_j\) se calcula como:
\[ IC_{95\%}(\beta_j) = \hat{\beta}_j \pm t^* \cdot SE(\hat{\beta}_j) \]
donde: - \(\hat{\beta}_j\) es la estimación puntual del coeficiente, - \(SE(\hat{\beta}_j)\) es el error estándar de dicha estimación, - \(t^*\) es el valor crítico de la distribución t de Student con \(n - p - 1\) grados de libertad (número de observaciones menos número de parámetros).
Un intervalo que incluye el cero sugiere que no hay evidencia suficiente para afirmar que el coeficiente sea distinto de cero al nivel de confianza elegido.
## 2.5 % 97.5 %
## (Intercept) 3.086557749 7.18841554
## cajas 1.245711004 1.80645228
## distancia 0.002363293 0.01465241
Es importante destacar que, aunque las observaciones \((X_{i1}, X_{i2}, \dots, X_{i(p-1)}, Y_i)\) no pueden representarse gráficamente en más de dos dimensiones, siempre es posible analizar gráficamente:
A continuación, se listan gráficos útiles para evaluar la idoneidad del modelo de regresión múltiple:
Recordemos que el leverage mide qué tan lejos están las observaciones respecto del centro de los datos.
\[h_{ii}>2\bar{h}\]
En regresión lineal simple, \(\sum_{i=1}^n h_{ii}=2\). Por lo tanto, \(\bar{h}=\dfrac{2}{n}.\) Entonces, una observación se considerará que tiene alto leverage si \[h_{ii}>\dfrac{4}{n}.\]
En regresión lineal múltiple, \(\sum_{i=1}^n h_{ii}=p\). Por lo tanto, \(\bar{h}=\dfrac{p}{n}.\) Entonces, una observación se considerará que tiene alto leverage si \[h_{ii}>\dfrac{2 p}{n}\]
Cabe señalar, que este criterio está ajustado por la cantidad de covariables.
Algunos autores utilizan otro criterio sin considerar este ajuste.
Se considera que una observación tiene alto leverage si \[h_{ii}>0.5.\]
Se considera que una observación tiene un moderado leverage si \[0.2<h_{ii}<0.5.\]
También se puede graficar un histograma de los valores \(h_{ii}\) y observar si hay algún valor/es que se comporta/n encuentran lejanos al resto.
## 1 2 3 4 5 6 7
## 0.11992977 0.10097319 0.10024969 0.10128434 0.16642759 0.08066357 0.07620000
## 8 9 10 11 12 13 14
## 0.06254612 0.06537738 0.11367920 0.07741039 0.06437975 0.04290146 0.09998709
## 15 16 17 18 19 20 21
## 0.04661503 0.04687996 0.08033747 0.11083391 0.20438000 0.13894064 0.14675966
## 22 23 24
## 0.18537865 0.21115081 0.55671434
## 24
## 24
Mide la influencia que tiene una observación determinada en la estimación de los parámetros de la regresión.
Una observación con distancia de Cook grande indica que la ausencia de esa observación podría cambiar los valores estimados de los parámetros.
Es una medida que toma en cuenta tanto los residuos como los puntos de alto leverage.
\[ D_i = \frac{\displaystyle{\sum_{j=1}^n}\left(\widehat{Y}_j - \widehat{Y}_{(i)j}\right)^2}{p \widehat{\sigma}^2} \] donde: - \(p\) es la cantidad de parámetros a estimar en el modelo.
\(\widehat{Y}_{(i)j}\) corresponde al valor predicho para la j-ésima observación si se usaron las \(n-1\) restantes observaciones para hacer el ajuste.
\(\widehat{Y}_{j}\) es el valor predicho para la j-ésima observación en el modelo ajustado con las \(n\) observaciones.
\[ D_i = \dfrac{r_{si} ^2}{p} \frac{h_{ii}}{1 - h_{ii}} \] donde \(r_{si}\) es el i-ésimo residuo standarizado y \(h_{ii}\) es el i-ésimo valor del leverage, ambos definidos anteriormente.
El \(r_{si}\) mide cuán lejos está el valor predicho de la observación, mientras que el \(h_{ii}\) mide el grado de apalancamiento de la observación \(i\). Así, un valor elevado de \(D_{i}\) puede deberse a un gran valor de \(r_{i}\) , a un gran valor de \(h_{ii}\) o a ambos.
Entonces:
\(D_i\) depende de dos factores:
los \(r_i\) y los \(h_{ii}\).
Cuánto más grande son los \(r_i\) y los \(h_{ii}\), mayor será el valor de \(D_i\).
Por lo tanto, la \(j-ésima\) observación puede ser influyente si sucede alguna de las siguientes situaciones:
tener un \(r_i\) alto y sólo un moderado valor de \(h_{ii}\).
tener un valor alto valor de \(h_{ii}\) con sólo un moderado valor de \(r_i\).
tener tanto un alto valor de \(r_i\) como un alto valor de \(r_i\).
Los valores de \(D_i\) se comparan con los percentiles de la distribución \(F\) de Fisher con \(p\) y \(n-p\) grados de libertad en el numerador y denominador, respectivamente.
El criterio para decidir si una observación es influyente es el siguiente:
Si \(D_i<\) percentil 0,20 de la distribución \(F_{p, n-p}\) entonces la observación no es influyente.
Si \(D_i>\) percentil 0,50 de la distribución \(F_{p, n-p}\) entonces la observación es muy influyente y requerirá tomar alguna medida.
Si \(D_i\) se encuentra entre el percentil 0,20 y el percentil 0,50 de la distribución \(F_{p, n-p}\) se sugiere mirar además otros estadisticos.
En el ejemplo:
## [1] 0.8148716
## [1] 0.3352225
## 1 2 3 4 5 6
## 2.221728e-03 5.753886e-02 8.259894e-03 2.793573e-02 1.716661e-01 1.664864e-02
## 7 8 9 10 11 12
## 1.469579e-02 7.105612e-02 1.964887e-04 2.393920e-02 3.662926e-05 3.045623e-05
## 13 14 15 16 17 18
## 6.607918e-04 2.174400e-03 4.678501e-03 4.241174e-02 1.163750e-02 1.125283e-01
## 19 20 21 22 23 24
## 2.184760e-01 8.405315e-02 9.061495e-02 9.429366e-02 7.143746e-02 5.606107e-02
Hay varias formas.
# Distancia de Cook
cooksd <- cooks.distance(modelo)
# Ver los índices de los puntos más influyentes
which(cooksd >limite05.Cook)## named integer(0)
## named integer(0)
##
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos.sin)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9335 -1.8143 -0.2151 2.0083 4.3673
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.347585 1.165595 4.588 0.000178 ***
## cajas 1.480844 0.186832 7.926 1.35e-07 ***
## distancia 0.008750 0.003093 2.829 0.010364 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.57 on 20 degrees of freedom
## Multiple R-squared: 0.8951, Adjusted R-squared: 0.8846
## F-statistic: 85.35 on 2 and 20 DF, p-value: 1.61e-10
ggplot(data = Tiempos.sin, aes(x = fitted(modelo.sin), y = rstandard(modelo.sin))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "valores ajustados", y = "residuos estandarizados",
title = "Residuos estandarizados vs Valores predichos")ggplot(data = Tiempos.sin, aes(x = cajas, y = rstandard(modelo.sin))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "cantidad de cajas", y = "residuos estandarizados",
title = "Residuos estandarizados vs Cantidad de cajas")ggplot(data = Tiempos.sin, aes(x = distancia, y = rstandard(modelo.sin))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "distancia", y = "residuos estandarizados",
title = "Residuos estandarizados vs Distancia")Normalidad de los residuos
residuos.sin.df=data.frame(resid=rstandard(modelo.sin))
ggplot(data = residuos.sin.df, aes(sample = rstandard(modelo.sin))) +
stat_qq() +
stat_qq_line(color = "red") +
labs(title = "QQ Plot de los residuos",
x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
theme_minimal()ggplot(residuos.sin.df, aes(y = rstandard(modelo.sin))) +
geom_boxplot(fill = "skyblue") +
labs(title = "Boxplot de residuos estandarizados",
y = "Residuos estandarizados") +
theme_minimal()No se detectan cambios grande en la estimación de los parámetros, la estimación de la varianza de los residuos, en las decisiones de los test individuales ni en el test global. Por lo tanto se elige el modelo con la observación \(24\).
library(ggplot2)
library(GGally)
ggpairs(Tiempos2,
columns = 1:3,
mapping = aes(alpha = 0.9),
upper = list(continuous = wrap("points", color = "darkred")),
lower = list(continuous = wrap("points", color = "blue")),
diag = list(continuous = wrap("densityDiag", fill = "lightblue")))##
## Call:
## lm(formula = tiempo ~ cajas + distancia, data = Tiempos2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7880 -0.6629 0.4364 1.1566 7.4197
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.341231 1.096730 2.135 0.044170 *
## cajas 1.615907 0.170735 9.464 3.25e-09 ***
## distancia 0.014385 0.003613 3.981 0.000631 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.259 on 22 degrees of freedom
## Multiple R-squared: 0.9596, Adjusted R-squared: 0.9559
## F-statistic: 261.2 on 2 and 22 DF, p-value: 4.687e-16
# Se crea el grafico 3d y se guarda en un objeto, por ejemplo mi_3d
mi_3d <- scatterplot3d(x=Tiempos2$cajas, y=Tiempos2$distancia, z=Tiempos2$tiempo, pch=16, cex.lab=1,highlight.3d=TRUE, type="h", xlab="Cantidad de cajas",ylab="Distancia ", zlab="Tiempo")
# Para agregar el plano usamos $plane3d( ) con argumento modelo ajustado
mi_3d$plane3d(modelo2, lty.box="solid", col="mediumblue")ggplot(data = Tiempos2, aes(x = fitted(modelo2), y = rstandard(modelo2))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "valores ajustados", y = "residuos estandarizados",
title = "Residuos estandarizados vs Valores predichos")Estos gráficos sirven para detectar la existencia de una relación curvilínea entre los residuos y cada variable regresora. Esto indicaría que debiera incluirse un término \(X_i^2\) o \(ln(X_i)\) o alguna otra transformación.
ggplot(data = Tiempos2, aes(x = cajas, y = rstandard(modelo2))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "cantidad de cajas", y = "residuos estandarizados",
title = "Residuos estandarizados vs Cantidad de cajas")ggplot(data = Tiempos2, aes(x = distancia, y = rstandard(modelo2))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "distancia", y = "residuos estandarizados",
title = "Residuos estandarizados vs Distancia")residuos.df=data.frame(resid=rstandard(modelo2))
ggplot(data = residuos.df, aes(sample = rstandard(modelo2))) +
stat_qq() +
stat_qq_line(color = "red") +
labs(title = "QQ Plot de los residuos",
x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
theme_minimal()## 20 1 24 22 23 21
## -1.87354643 -1.62767993 -1.49605875 -1.44999541 -1.44368977 -0.87784258
## 16 12 5 6 25 3
## -0.22270023 -0.19325733 -0.14176094 -0.09080847 -0.06750861 -0.01609165
## 17 15 7 13 14 2
## 0.13803929 0.21029137 0.27042496 0.32517935 0.34113547 0.36484267
## 8 19 11 10 18 4
## 0.36672118 0.57876634 0.71807970 0.81325432 1.11295196 1.57972040
## 9
## 3.21376278
## 1 2 3 4 5 6 7
## 0.11992977 0.10097319 0.10024969 0.10128434 0.16642759 0.08066357 0.07620000
## 8 9 10 11 12 13 14
## 0.06254612 0.06537738 0.11367920 0.07741039 0.06437975 0.04290146 0.09998709
## 15 16 17 18 19 20 21
## 0.04661503 0.04687996 0.08033747 0.11083391 0.20438000 0.13894064 0.14675966
## 22 23 24
## 0.18537865 0.21115081 0.55671434
## 24
## 0.5567143
## named numeric(0)
La multicolinealidad ocurre cuando existe correlación entre dos o más variables regresoras.
Relaciones lógicas o matemáticas entre
variables
Variables que están construidas a partir de otras.
Ejemplo: población total y población urbana.
Variables derivadas unas de otras
Incluir tanto una variable como una versión transformada o
combinada.
Ejemplo: incluir tanto peso como el índice de masa corporal
(IMC).
Variables naturalmente correlacionadas por
contexto
Algunas variables tienden a moverse juntas por razones sociales o
económicas.
Ejemplo: nivel educativo y salario.
Redundancia por variables similares
Incluir indicadores que miden aspectos muy parecidos.
Ejemplo: gasto mensual y gasto anual.
Tamaño de muestra pequeño
Con pocos datos, las correlaciones entre variables explicativas se
amplifican, generando efectos de colinealidad incluso con relaciones
moderadas.
Los estimadores de los coeficientes pueden cambiar drásticamente con pequeños cambios en los datos.
Cuando se incluyen covariables altamente correlacionadas en un modelo, los errores estándar de los coeficientes estimados tienden a aumentar de manera artificial. A este fenómeno se lo conoce como inflación de la varianza de los estimadores.
Los coeficientes pueden no resultar estadísticamente significativos, incluso cuando existe una relación real entre la variable dependiente y las variables explicativas.
Un método para detectar colinealidad entre las variables regresoras es el Factor de Infalción de la Varianza (VIF por sus siglas en inglés).
Este factor se calcula para cada variable regresora.
El VIF para la k-ésima variable regresora se define:
\[\mathrm{VIF}_k=\dfrac{1}{1-R_k^2} \] con \(k=1,\ldots,p-1\) y \(R_k^2\) es el coeficiente de determinación múltiple cuando \(X_k\) es regresado en las \(p- 2\) restantes covariables \(X\) en el modelo.
\(\mathrm{VIF}_k=1\) cuando \(R_k^2=0\) : la covariable \(X_k\) no está correlacionada con las demás.
Si \(R_k^2>0\), entonces \(\mathrm{VIF}_k>1\) : hay algún grado de colinealidad.
Si \(R_k^2 \approx 1\), el \(\mathrm{VIF}_k\) puede ser muy grande, indicando alta colinealidad.
\(\mathrm{VIF} > \mathbf{10} \rightarrow\) indicio fuerte de multicolinealidad.
\(\mathbf{5}<\mathrm{VIF} < \mathbf{10}\) indicio de moderada multicolinealidad.
Otro criterio: si el promedio de los VIF es significativamente mayor a 1, también sugiere problemas de colinealidad.
A continuación se muestra un ejemplo comparando dos modelos de regresión:
Uno sin multicolinealidad, con tres variables independientes no correlacionadas.
Otro con multicolinealidad, donde una variable es casi una copia de otra.
set.seed(123)
# N° de observaciones
n <- 100
# Modelo sin multicolinealidad
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y1 <- 3 + 2 * x1 + 1.5 * x2 + 1 * x3 + rnorm(n)
datos1 <- data.frame(y = y1, x1 = x1, x2 = x2, x3 = x3)
# Modelo con multicolinealidad (x3 ≈ x1)
x3c <- x1 + rnorm(n, 0, 0.05) # Alta correlación con x1
y2 <- 3 + 2 * x1 + 1.5 * x2 + 1 * x3c + rnorm(n)
y2 <- 3 + 2 * x1 + 1.5 * x2 + 1 * x3c + rnorm(n)
datos2 <- data.frame(y = y2, x1 = x1, x2 = x2, x3 = x3c)##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = datos1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.49138 -0.65392 0.05664 0.67033 2.53210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.9807 0.1073 27.768 < 2e-16 ***
## x1 1.9445 0.1169 16.638 < 2e-16 ***
## x2 1.5462 0.1095 14.126 < 2e-16 ***
## x3 0.9426 0.1122 8.399 4.04e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.052 on 96 degrees of freedom
## Multiple R-squared: 0.8392, Adjusted R-squared: 0.8342
## F-statistic: 167 on 3 and 96 DF, p-value: < 2.2e-16
ggplot(data = datos1 , aes(x = fitted(modelo_sin_col), y = rstandard(modelo_sin_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "valores ajustados", y = "residuos estandarizados",
title = "Residuos estandarizados vs Valores predichos")ggplot(data = datos1 , aes(x = x1, y = rstandard(modelo_sin_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = expression(x[1]), y = "residuos estandarizados",
title = expression(paste("Residuos estandarizados vs. ", x[1])))ggplot(data = datos1 , aes(x = x2, y = rstandard(modelo_sin_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = expression(x[2]), y = "residuos estandarizados",
title = expression(paste("Residuos estandarizados vs. ", x[2])))ggplot(data = datos1 , aes(x = x3, y = rstandard(modelo_sin_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = expression(x[3]), y = "residuos estandarizados",
title = expression(paste("Residuos estandarizados vs. ", x[3])))residuos.df=data.frame(resid=rstandard(modelo_sin_col))
ggplot(data = residuos.df, aes(sample = rstandard(modelo_sin_col))) +
stat_qq() +
stat_qq_line(color = "red") +
labs(title = "QQ Plot de los residuos",
x = "Cuantiles teóricos", y = "Cuantiles de los residuos") +
theme_minimal()##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = datos2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.43303 -0.74973 0.04546 0.70439 2.40904
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8625 0.1058 27.047 <2e-16 ***
## x1 2.0176 2.1563 0.936 0.352
## x2 1.5811 0.1094 14.452 <2e-16 ***
## x3 0.9474 2.1768 0.435 0.664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.04 on 96 degrees of freedom
## Multiple R-squared: 0.8974, Adjusted R-squared: 0.8942
## F-statistic: 280 on 3 and 96 DF, p-value: < 2.2e-16
ggplot(data = datos2, aes(x = fitted(modelo_con_col), y = rstandard(modelo_con_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = "valores ajustados", y = "residuos estandarizados",
title = "Residuos estandarizados vs Valores predichos")ggplot(data = datos2, aes(x = x1, y = rstandard(modelo_con_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = expression(x[1]), y = "residuos estandarizados",
title = expression(paste("Residuos estandarizados vs. ", x[1])))ggplot(data = datos2, aes(x = x2, y = rstandard(modelo_con_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = expression(x[2]), y = "residuos estandarizados",
title = expression(paste("Residuos estandarizados vs. ", x[2])))ggplot(data = datos2, aes(x = x3, y = rstandard(modelo_con_col))) +
geom_point() +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(x = expression(x[3]), y = "residuos estandarizados",
title = expression(paste("Residuos estandarizados vs. ", x[3])))residuos.df=data.frame(resid=rstandard(modelo_con_col))
ggplot(data = residuos.df, aes(sample = rstandard(modelo_con_col))) +
stat_qq() +
stat_qq_line(color = "red") +
labs(title = "QQ Plot de los residuos",
x = "cuantiles teóricos", y = "cuantiles de los residuos") +
theme_minimal()## x1 x2 x3
## 1.019125 1.003057 1.017576
## x1 x2 x3
## 354.260148 1.023393 354.548013
Lo más fácil es elegir un subconjunto de las variables regresoras poco correlacionadas.
El asunto es… si detectamos dos variables muy correlacionadas ¿cuál omito? En general, conviene omitir aquella variable que tenga:
Otra posibilidad es utilizar métodos para eliminar variables a través de procedimientos de selección automáticos (se presentarán más adelante).
También se puede consturir una tercer variable que sea combinación de estas variables correlacionadas. Por ejemplo: el índice de precios.
Selección de variables: con este método trabajaremos.
Métodos de regularización o penalización. Ridge o Lasso: mejoran la estabilidad de los coeficientes, especialmente cuando hay multicolinealidad o muchas variables. Añaden un término de penalización al criterio de ajuste del modelo. Es decir, en lugar de minimizar solo el error cuadrático (como en la regresión lineal clásica), también se penalizan los coeficientes grandes.
Métodos de reducción de la dimensión: son técnicas que reducen la cantidad de variables regresoras proyectando sobre un espacio de dimensión menor. Este método no se tratará en el curso.
Reformulamos el modelo eligiendo \(x_1\) o \(x_3\).
Referencias:
Montgomery, E., Vining, D. y Peck. (2006). Introducción Al Análisis de Regresión Lineal. 3ed ed. México: Cecsa.
James, G., Witten, D., Hastie, T., y Tibshirani, R. (2013). An introduction to statistical learning (Vol. 112). Springer.